2. Descriptive Statistics

Means SDs Skewness Kurtosis Betas

Mean (price) Mean(return %) SD (price) SD(return %) Skewness (price) Skewness (return %) Kurtosis (price) Kurtosis (return %) Beta
MSFT 117.508 2.495 86.817 5.918 1.043 0.231 -0.080 0.691 1.004
TSLA 203.202 4.802 298.245 17.296 1.823 1.293 1.974 2.342 1.993
AAPL 58.870 2.774 45.183 7.744 1.271 -0.233 0.292 -0.215 1.247
TWTR 34.921 0.506 14.229 14.535 0.650 0.416 0.057 0.446 0.869
AMZN 1550.228 2.595 1062.254 8.163 0.549 0.402 -1.051 0.727 1.175
FB 171.007 1.636 80.593 8.123 0.713 -0.335 -0.240 2.892 1.219
NFLX 264.865 2.569 176.639 11.741 0.454 0.472 -0.997 1.115 0.941
AAL 33.150 0.112 11.580 11.493 -0.356 -0.096 -0.899 0.284 1.559
DAL 42.931 0.840 8.403 9.510 -0.057 -0.187 -0.600 2.746 1.183
BAC 23.756 1.398 9.369 8.191 0.696 -0.159 -0.265 1.109 1.438
NVDA 64.032 5.146 74.026 11.778 1.699 0.071 2.263 0.489 1.452
WBD 29.026 0.106 6.327 11.332 1.093 0.854 1.452 1.192 1.244
S.P500 2770.638 1.036 799.620 4.002 0.966 -0.393 -0.052 1.433 1.000

Plots

Equity curve

Stationary Test

MSFT 0.1
TSLA 0.1
AAPL 0.1
TWTR 0.1
AMZN 0.1
FB 0.1
NFLX 0.1
AAL 0.1
DAL 0.1
BAC 0.1
NVDA 0.1
WBD 0.1
S&P500 0.1

Hist, Boxplot, qqplot

Distributions

Asset t-distribution Normal Distribution GED
MSFT -139.66549 -271.25573 -276.23424
TSLA -38.58013 -80.34437 -69.20914
AAPL -112.15671 -220.30961 -218.90939
TWTR -50.63192 -97.29781 -94.90164
AMZN -107.95913 -209.86197 -209.57404
FB -111.98310 -209.20876 -217.11195
NFLX -73.72148 -138.40822 -140.89449
AAL -73.60630 -140.92136 -141.25665
DAL -96.76172 -178.09420 -184.85237
BAC -108.18352 -207.92323 -209.41384
NVDA -71.69179 -136.35263 -137.27907
WBD -77.72453 -151.06318 -153.27377

Sharpe’s Slope

Assets Sharpe’s Slope
MSFT 0.3012400
TSLA 0.2364437
AAPL 0.2662491
TWTR -0.0141628
AMZN 0.2307159
FB 0.1137256
NFLX 0.1581712
AAL -0.0522252
DAL 0.0134156
BAC 0.0838003
NVDA 0.3764550
WBD -0.0534724

M to Y

Monthly Mean Monthly SD Annual Mean Annual SD
MSFT 2.4946532 5.917747 29.935838 20.49968
TSLA 4.8015901 17.296293 57.619081 59.91612
AAPL 2.7737763 7.743821 33.285316 26.82538
TWTR 0.5061345 14.535005 6.073614 50.35073
AMZN 2.5953512 8.163112 31.144214 28.27785
FB 1.6357391 8.122604 19.628869 28.13753
NFLX 2.5691372 11.741366 30.829647 40.67329
AAL 0.1117527 11.493264 1.341033 39.81384
DAL 0.8395770 9.510284 10.074924 32.94459
BAC 1.3984216 8.191269 16.781059 28.37539
NVDA 5.1459413 11.778168 61.751296 40.80077
WBD 0.1060357 11.332111 1.272428 39.25558
S&P500 1.0358383 4.002103 12.430059 13.86369

Pairewise

Covariance Matrix

MSFT TSLA AAPL TWTR AMZN FB NFLX AAL DAL BAC NVDA WBD S&P500
MSFT 0.0035 0.0038 0.0024 0.0005 0.0027 0.0019 0.0026 0.0020 0.0014 0.0021 0.0034 0.0020 0.0016
TSLA 0.0038 0.0299 0.0063 0.0036 0.0045 0.0041 0.0055 0.0032 0.0020 0.0034 0.0059 0.0043 0.0032
AAPL 0.0024 0.0063 0.0060 0.0024 0.0028 0.0030 0.0025 0.0024 0.0019 0.0016 0.0047 0.0019 0.0020
TWTR 0.0005 0.0036 0.0024 0.0211 0.0023 0.0029 0.0018 0.0012 0.0023 0.0024 0.0019 0.0050 0.0014
AMZN 0.0027 0.0045 0.0028 0.0023 0.0067 0.0032 0.0056 0.0015 0.0008 0.0020 0.0042 0.0017 0.0019
FB 0.0019 0.0041 0.0030 0.0029 0.0032 0.0066 0.0036 0.0019 0.0008 0.0023 0.0028 0.0021 0.0019
NFLX 0.0026 0.0055 0.0025 0.0018 0.0056 0.0036 0.0138 0.0013 0.0001 0.0018 0.0051 0.0021 0.0015
AAL 0.0020 0.0032 0.0024 0.0012 0.0015 0.0019 0.0013 0.0132 0.0086 0.0050 0.0032 0.0039 0.0025
DAL 0.0014 0.0020 0.0019 0.0023 0.0008 0.0008 0.0001 0.0086 0.0090 0.0039 0.0024 0.0036 0.0019
BAC 0.0021 0.0034 0.0016 0.0024 0.0020 0.0023 0.0018 0.0050 0.0039 0.0067 0.0027 0.0043 0.0023
NVDA 0.0034 0.0059 0.0047 0.0019 0.0042 0.0028 0.0051 0.0032 0.0024 0.0027 0.0139 0.0008 0.0023
WBD 0.0020 0.0043 0.0019 0.0050 0.0017 0.0021 0.0021 0.0039 0.0036 0.0043 0.0008 0.0128 0.0020
S&P500 0.0016 0.0032 0.0020 0.0014 0.0019 0.0019 0.0015 0.0025 0.0019 0.0023 0.0023 0.0020 0.0016

3. Portfolio Theory

With Short Sale

library(quadprog)
R = 100*return[,2:13]
mean_p <- apply(R,2,mean)
cov_p <- cov(R)
sd_vect_p <- sqrt(diag(cov_p))
# min(mean_p)
# max(mean_p)
### With shortsale
M_p = length(mean_p)
Amat_p <- cbind(rep(1,M_p),mean_p)
mu_P = seq(0.07, 5.4, length = 300)
# Target portfolio means for the expect portfolio return
sd_P = mu_P # set up storage for std dev's of portfolio returns
weights_p = matrix(0, nrow = 300, ncol = M_p) # storage for return
for (i in 1:length(mu_P)) { # find the optimal portfolios
  bvec_p <- c(1, mu_P[i])
  result_p = solve.QP(Dmat = 2 * cov_p, dvec = rep(0, M_p), Amat = Amat_p, 
                    bvec = bvec_p, meq = 2)
  sd_P[i] = sqrt(result_p$value)
  weights_p[i, ] = result_p$solution
}
plot(sd_P, mu_P, type = "l", xlim = c(0,15), ylim = c(0, 6), lty = 3, lwd = 2, main = "With Short Sale") 
# plot efficient frontier (and inefficient portfolios below the min var portfolio)
mufree_p = mean(return$`Treasury Bill 3 month (rf)`)# input value of risk-free interest rate
points(0, mufree_p, cex = 4, pch = "*") # show risk-free asset
sharpe_p = (mu_P - mufree_p) / sd_P # compute Sharpes ratios
ind_p = (sharpe_p == max(sharpe_p)) # Find maximum Sharpes ratio
#weights_p[ind_p,] # print the weights of the tangency portfolio
lines(c(0, 15), mufree_p + c(0, 15) * (mu_P[ind_p] - mufree_p) / sd_P[ind_p], lwd = 4, 
      lty = 1, col = "blue") # show line of optimal portfolios
points(sd_P[ind_p], mu_P[ind_p], cex = 4, pch = "*") # tangency portfolio
ind2_p = (sd_P == min(sd_P)) # find the minimum variance portfolio
points(sd_P[ind2_p], mu_P[ind2_p], cex = 2, pch = "+") # min var portfolio
ind3_p = (mu_P > mu_P[ind2_p])
lines(sd_P[ind3_p], mu_P[ind3_p], type = "l", xlim = c(0, 25), ylim = c(0,30),
      lwd = 3, col = 'red') # plot the efficient frontier
for(i in 1:12){
  text(sd_vect_p[i], mean_p[i],colnames(return[,i+1]), cex=0.8)
}

### MVP
(mvp_meanreturn <- mu_P[ind2_p])
## [1] 1.816957
(mvp_sd <- sd_P[ind2_p])
## [1] 4.980327
weights_mvp <- weights_p[ind2_p,]
weights_mvp <- t(data.frame(weights_mvp))
colnames(weights_mvp) <- colnames(return[2:13])
knitr::kable(round(weights_mvp*100,2))
MSFT TSLA AAPL TWTR AMZN FB NFLX AAL DAL BAC NVDA WBD
weights_mvp 50.17 -4.98 13.48 4.89 3.85 13.16 4.31 -9.16 21.82 8.07 -5.96 0.34
(mvp_meanreturn_ann <- mvp_meanreturn*12)
## [1] 21.80348
(mvp_sd_ann <- mvp_sd*sqrt(12))
## [1] 17.25236
### Efficient Portfolio Frontier
EPF_mean <- mu_P[ind3_p]
EPF_sd <- sd_P[ind3_p]

### Tangency Portfolio
(tan_meanreturn <- mu_P[ind_p])
## [1] 5.4
(tan_sd <- sd_P[ind_p])
## [1] 9.86121
(tan_var <- tan_sd^2)
## [1] 97.24347
(tan_sharpes <- (tan_meanreturn-mufree_p)/tan_sd)
## [1] 0.4753989

VaR&ES

### Tail dependence can be seen among the assets, therefore we can fit
### our portfolio with multivariate t-distribution.

## MVP VaR&ES
library(MASS)
alpha = 0.05
return_mvp <- rowSums(data.frame(
  weights_mvp[1] * return[, 2],
  weights_mvp[2] * return[, 3],
  weights_mvp[3] * return[, 4],
  weights_mvp[4] * return[, 5],
  weights_mvp[5] * return[, 6],
  weights_mvp[6] * return[, 7],
  weights_mvp[7] * return[, 8],
  weights_mvp[8] * return[, 9],
  weights_mvp[9] * return[, 10],
  weights_mvp[10] * return[, 11],
  weights_mvp[11] * return[, 12],
  weights_mvp[12] * return[, 13]))
fitt_mvp = fitdistr(return_mvp,"t")
param_mvp = as.numeric(fitt_mvp$estimate)
mean_mvpfit = param_mvp[1]
df_mvpfit = param_mvp[3]
sd_mvpfit = param_mvp[2] * sqrt((df_mvpfit) / (df_mvpfit - 2))
lambda_mvpfit = param_mvp[2]
qalpha_mvp = qt(alpha, df = df_mvpfit)
VaR_par_mvp = -100000 * (mean_mvpfit + lambda_mvpfit * qalpha_mvp)
es1_mvp = dt(qalpha_mvp, df = df_mvpfit) / (alpha)
es2_mvp=(df_mvpfit+qalpha_mvp^2)/(df_mvpfit-1)
es3_mvp=-mean_mvpfit+lambda_mvpfit*es1_mvp*es2_mvp
ES_par_mvp = 100000*es3_mvp
VaR_par_mvp
## [1] 6285.691
ES_par_mvp
## [1] 8957.941
## Asset VaR
S0 = 100000
qnalpha = qnorm(0.05)

### MSFT
q_msft = as.numeric(quantile(return$MSFT, alpha))
VAR_msft = -S0 * q_msft
#VAR_msft

### TSLA
fit_tsla <- fitdistr(return$TSLA, "normal")
param_tsla = as.numeric(fit_tsla$estimate)
mean_tsla = param_tsla[1]
sd_tsla = param_tsla[2]
VAR_tsla = -S0*(mean_tsla+qnalpha*sd_tsla)
#VAR_tsla

### AAPL
fit_aapl <- fitdistr(return$AAPL, "normal")
param_aapl = as.numeric(fit_aapl$estimate)
mean_aapl = param_aapl[1]
sd_aapl = param_aapl[2]
VAR_aapl = -S0*(mean_aapl+qnalpha*sd_aapl)
#VAR_aapl

### TWTR
fit_twtr <- fitdistr(return$TWTR, "normal")
param_twtr = as.numeric(fit_twtr$estimate)
mean_twtr = param_twtr[1]
sd_twtr = param_twtr[2]
VAR_twtr = -S0*(mean_twtr+qnalpha*sd_twtr)
#VAR_twtr

### AMZN
fit_amzn <- fitdistr(return$AMZN, "normal")
param_amzn = as.numeric(fit_amzn$estimate)
mean_amzn = param_amzn[1]
sd_amzn = param_amzn[2]
VAR_amzn = -S0*(mean_amzn+qnalpha*sd_amzn)
#VAR_amzn

### FB
q_fb = as.numeric(quantile(return$FB, alpha))
VAR_fb = -S0 * q_fb
#VAR_fb

### NFLX
q_nflx = as.numeric(quantile(return$NFLX, alpha))
VAR_nflx = -S0 * q_nflx
#VAR_nflx

### AAL
q_aal = as.numeric(quantile(return$AAL, alpha))
VAR_aal = -S0 * q_aal
#VAR_aal

### DAL
q_dal = as.numeric(quantile(return$DAL, alpha))
VAR_dal = -S0 * q_dal
#VAR_dal

### BAC
q_bac = as.numeric(quantile(return$BAC, alpha))
VAR_bac = -S0 * q_bac
#VAR_bac

### NVDA
q_nvda = as.numeric(quantile(return$NVDA, alpha))
VAR_nvda = -S0 * q_nvda
#VAR_nvda

### WBD
q_wbd = as.numeric(quantile(return$WBD, alpha))
VAR_wbd = -S0 * q_wbd
#VAR_wbd

VAR_asset <- c(VAR_msft, VAR_tsla, VAR_aapl, VAR_twtr, VAR_amzn, VAR_fb, VAR_nflx,
               VAR_aal, VAR_dal, VAR_bac, VAR_nvda, VAR_wbd)
var_df <- cbind(names, VAR_asset)
var_df <- rbind(var_df,VaR_par_mvp)
var_df[13,1] <- c("MVP")
var_df[13,2] <- 6285.69101812641
colnames(var_df) <- c("Assets", "Var")

knitr::kable(var_df)
Assets Var
MSFT 6918.17869693084
TSLA 23502.7562034031
AAPL 9898.52221650451
TWTR 23279.5289061133
AMZN 10763.0915501661
FB 10561.0446590248
NFLX 13072.710950181
AAL 15706.0583956349
DAL 12891.9240002285
BAC 13041.1530377474
NVDA 11974.7969532429
WBD 14903.9099722124
VaR_par_mvp MVP 6285.69101812641

Assets’ Sharpe’s Ratios

names_sh <- data.frame(colnames(return[1,2:14]))
sharpes_list <- rep(NA, 13)
for(i in 2:14){
  sharpes_list[i-1] = (mean(unlist(return[,i]))-mean(unlist(return[,15]))/100)/sd(unlist(return[,i]))
}
sharpes_list <- data.frame(sharpes_list)
shar_df <- cbind(names_sh,sharpes_list)
shar_df
##    colnames.return.1..2.14.. sharpes_list
## 1                       MSFT   0.30124000
## 2                       TSLA   0.23644366
## 3                       AAPL   0.26624907
## 4                       TWTR  -0.01416282
## 5                       AMZN   0.23071594
## 6                         FB   0.11372558
## 7                       NFLX   0.15817121
## 8                        AAL  -0.05222523
## 9                        DAL   0.01341556
## 10                       BAC   0.08380025
## 11                      NVDA   0.37645499
## 12                       WBD  -0.05347243
## 13                    S&P500   0.08091924

Without shortshale

### Without shortsale
#R = 100*return[,2:13]
#mean_p <- apply(R,2,mean)
#cov_p <- cov(R)
#sd_vect_p <- sqrt(diag(cov_p))
### With shortsale
#M_p = length(mean_p)
Amat_p_noss <- cbind(rep(1,M_p),mean_p, diag(1,nrow=M_p))
mu_P_noss = seq(min(mean_p)+0.0001, max(mean_p)-0.0001, length = 300)
# Target portfolio means for the expect portfolio return
sd_P_noss = mu_P_noss # set up storage for std dev's of portfolio returns
weights_p_noss = matrix(0, nrow = 300, ncol = M_p) # storage for return
for (i in 1:length(mu_P_noss)) { # find the optimal portfolios
  bvec_p_noss <- c(1, mu_P_noss[i], rep(0,M_p))
  result_noss = solve.QP(Dmat = 2 * cov_p, dvec = rep(0, M_p), Amat = Amat_p_noss, 
                    bvec = bvec_p_noss, meq = 2)
  sd_P_noss[i] = sqrt(result_noss$value)
  weights_p_noss[i, ] = result_noss$solution
}
plot(sd_P_noss, mu_P_noss, type = "l", lty = 3, 
     lwd = 2, xlim = c(0,15), ylim = c(0,7), main = "Without Short Sale") 
# plot efficient frontier (and inefficient portfolios below the min var portfolio)
#mufree_p = mean(return$`Treasury Bill 3 month (rf)`) # input value of risk-free interest rate
points(0, mufree_p, cex = 4, pch = "*") # show risk-free asset
sharpe_p_noss = (mu_P_noss - mufree_p) / sd_P_noss # compute Sharpes ratios
ind_p_noss = (sharpe_p_noss == max(sharpe_p_noss)) # Find maximum Sharpes ratio
#weights_p[ind_p,] # print the weights of the tangency portfolio
lines(c(0, 15), mufree_p + c(0, 15) * (mu_P_noss[ind_p_noss] - mufree_p) / sd_P_noss[ind_p_noss], lwd = 4, 
      lty = 1, col = "blue") # show line of optimal portfolios
points(sd_P_noss[ind_p_noss], mu_P_noss[ind_p_noss], cex = 4, pch = "*") # tangency portfolio
ind2_p_noss = (sd_P_noss == min(sd_P_noss)) # find the minimum variance portfolio
points(sd_P_noss[ind2_p_noss], mu_P_noss[ind2_p_noss], cex = 2, pch = "+") # min var portfolio
ind3_p_noss = (mu_P_noss > mu_P_noss[ind2_p_noss])
lines(sd_P_noss[ind3_p_noss], mu_P_noss[ind3_p_noss], type = "l",
      lwd = 3, col = 'red') # plot the efficient frontier
for(i in 1:12){
  text(sd_vect_p[i], mean_p[i],colnames(return[,i+1]), cex=0.8)
}

### MVP
(mvp_meanreturn_noss <- mu_P_noss[ind2_p_noss])
## [1] 1.977063
(mvp_sd_noss <- sd_P_noss[ind2_p_noss])
## [1] 5.110456
weights_mvp_noss <- weights_p_noss[ind2_p_noss,]
weights_mvp_noss <- t(data.frame(weights_mvp_noss))
colnames(weights_mvp_noss) <- colnames(return[2:13])
weights_mvp_noss
##                      MSFT          TSLA       AAPL       TWTR       AMZN
## weights_mvp_noss 0.487966 -1.764314e-18 0.06592605 0.05609903 0.03247623
##                         FB       NFLX           AAL      DAL        BAC
## weights_mvp_noss 0.1346882 0.02256666 -2.931383e-17 0.146045 0.04453321
##                          NVDA         WBD
## weights_mvp_noss 6.441445e-18 0.009699598
(mvp_meanreturn_ann_noss <- mvp_meanreturn_noss*12)
## [1] 23.72476
(mvp_sd_ann_noss <- mvp_sd_noss*sqrt(12))
## [1] 17.70314
sum(weights_mvp_noss)
## [1] 1
rownames(weights_mvp_noss) <- c("weight")
knitr::kable(round(weights_mvp_noss*100,4))
MSFT TSLA AAPL TWTR AMZN FB NFLX AAL DAL BAC NVDA WBD
weight 48.7966 0 6.5926 5.6099 3.2476 13.4688 2.2567 0 14.6045 4.4533 0 0.97
### Efficient Portfolio Frontier
EPF_mean_noss <- mu_P_noss[ind3_p_noss]
EPF_sd_noss <- sd_P_noss[ind3_p_noss]

### Tangency Portfolio
(tan_meanreturn_noss <- mu_P_noss[ind_p_noss])
## [1] 3.999688
(tan_sd_noss <- sd_P_noss[ind_p_noss])
## [1] 7.96967
(tan_var_noss <- tan_sd_noss^2)
## [1] 63.51564
(tan_sharpes_noss <- (tan_meanreturn_noss-mufree_p)/tan_sd_noss)
## [1] 0.412526
alpha = 0.05
return_mvp_noss <- rowSums(data.frame(
  weights_mvp_noss[1] * return[, 2],
  weights_mvp_noss[2] * return[, 3],
  weights_mvp_noss[3] * return[, 4],
  weights_mvp_noss[4] * return[, 5],
  weights_mvp_noss[5] * return[, 6],
  weights_mvp_noss[6] * return[, 7],
  weights_mvp_noss[7] * return[, 8],
  weights_mvp_noss[8] * return[, 9],
  weights_mvp_noss[9] * return[, 10],
  weights_mvp_noss[10] * return[, 11],
  weights_mvp_noss[11] * return[, 12],
  weights_mvp_noss[12] * return[, 13]))
fitt_mvp_noss = fitdistr(return_mvp_noss,"t")
param_mvp_noss = as.numeric(fitt_mvp_noss$estimate)
mean_mvpfit_noss = param_mvp_noss[1]
df_mvpfit_noss = param_mvp_noss[3]
sd_mvpfit_noss = param_mvp_noss[2] * sqrt((df_mvpfit_noss) / (df_mvpfit_noss - 2))
lambda_mvpfit_noss = param_mvp_noss[2]
qalpha_mvp_noss = qt(alpha, df = df_mvpfit_noss)
VaR_par_mvp_noss = -100000 * (mean_mvpfit_noss + lambda_mvpfit_noss * qalpha_mvp_noss)
es1_mvp_noss = dt(qalpha_mvp_noss, df = df_mvpfit_noss) / (alpha)
es2_mvp_noss=(df_mvpfit_noss+qalpha_mvp_noss^2)/(df_mvpfit_noss-1)
es3_mvp_noss=-mean_mvpfit_noss+lambda_mvpfit_noss*es1_mvp_noss*es2_mvp_noss
ES_par_mvp_noss = 100000*es3_mvp_noss
VaR_par_mvp_noss
## [1] 6291.464
ES_par_mvp_noss
## [1] 9034.092

##4. Asset Allocation

Efficient Portfolio

#efficient portfolio
new_dat = return[,2:13]
return_year = 12*colMeans(new_dat)

cov_mat = cov(new_dat)
eff_port = efficient.portfolio(er = return_year,cov_mat,target.return = 0.06,shorts = FALSE)
eff_port$er/12
## [1] 0.004999994
#check return
ret=0
for (i in 1:12)
{
  val = as.numeric(return_year[i])*as.numeric(eff_port$weights[i])
  ret = ret+val
  
}
ret
## [1] 0.05999993
#monthly risk
eff_risk = eff_port$sd/sqrt(12)
eff_risk
## [1] 0.02134934
# vaR = -S0*(u + qnorm(0.05)*sd)
S0 = 100000
q = qnorm(0.05)
eff_var = -S0*(eff_port$er/12 + q*eff_risk)
eff_var
## [1] 3011.654
#excel
cut_off = qnorm(0.05,mean = (1+eff_port$er/12)*S0,sd = S0*eff_risk)
var2 = S0-cut_off
var2
## [1] 3011.654
var2/S0
## [1] 0.03011654
#Expected Shortfall
S0*((-eff_port$er/12) + (dnorm(qnorm(0.05))/0.05)* eff_risk)
## [1] 3903.755

Tangency portfolio

risk_free = colMeans(return[,15]/100)

tan_port = tangency.portfolio(er=return_year, cov.mat=cov_mat, risk.free=risk_free,shorts=FALSE)
tan_port
## Call:
## tangency.portfolio(er = return_year, cov.mat = cov_mat, risk.free = risk_free, 
##     shorts = FALSE)
## 
## Portfolio expected return:     0.4214667 
## Portfolio standard deviation:  0.06885209 
## Portfolio weights:
##   MSFT   TSLA   AAPL   TWTR   AMZN     FB   NFLX    AAL    DAL    BAC   NVDA 
## 0.4830 0.0686 0.0903 0.0000 0.0452 0.0000 0.0000 0.0000 0.0000 0.0000 0.3129 
##    WBD 
## 0.0000
# w1 for tangency portfolio, w1 for risk free
# w1+w2 = 1
# 0.4799845 * w1 + 0.08543894 * w2 = 0.06

# 1.3041868% for the tangency portfolio and -0.3041868% for the risk free

A = matrix(c(1, tan_port$er, 1, as.numeric(risk_free)), ncol=2)
B = matrix(c(1, 0.06),ncol = 1)
solve(A)%*%B
##           [,1]
## [1,] 0.1276228
## [2,] 0.8723772
#vaR
S0 = 100000
q = qnorm(0.05)
eff_var = -S0*(0.05/12 + q*0.06885209/sqrt(12))
eff_var
## [1] 2852.626
eff_var/S0
## [1] 0.02852626
#ES 
S0*((-0.05/12) + (dnorm(qnorm(0.05))/0.05)* 0.06885209/sqrt(12))
## [1] 3683.158

5. PCA

correlation matrix of the returns on your assets.

AAL and DAL are the assets that most highly correlated with correlation = 0.78. If we only consider about the 12 assets, then NFLX and DAL has the lowest correlation = 0.0089. The NFLX also has a correlation with Treasury Bill 3 month (risk free) = 0.006.

Diversification will reduce risk for the assets that are negatively correlated. For example, if the correlation between 2 assets is -1, which is the perfect value for diversification, if one asset goes up 5%, the other asset will goes down 5%. So combining these 2 assets would minimize our risk to 0. In our 12 assets, we did not see any negative correlation. As a result, diversification will not reduce risk.

(note: here are a few assets that are slightly correlated such as NFLX and DAL, NFLX and AAL, TWTR and AAL, NVDA and WBD, MSFT and TWTR.)

library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:psych':
## 
##     logit
## The following object is masked from 'package:fBasics':
## 
##     densityPlot
df <- return[,2:13]
# correlation matrix 
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...) {
    usr <- par("usr")
    on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    Cor <- abs(cor(x, y)) 
    txt <- paste0(prefix, format(c(Cor, 0.123456789), digits = digits)[1])
    if(missing(cex.cor)) {
        cex.cor <- 0.4 / strwidth(txt)
    }
    text(0.5, 0.5, txt,
         cex = 1 + cex.cor * Cor) 
}

pairs(df,
      upper.panel = panel.cor,    # Correlation panel
      lower.panel = panel.smooth) # Smoothed regression lines

Parallel coordinate plot

parcoords(df, rownames = F, reorderable = T , queue = T, 
          withD3 = TRUE, alpha=0.5)

PCA Analysis

Another way to dig into the correlation between the return of the assets, we looked at PCA and its plots. PCA helps to reduce the dimension of the dataset and compute principal components that contain most of the information in the data. From the PCA Analysis, we found that first eight components captures almost 90% of the total variance, while nine components could explain 93% of the variance. The Scree Plot shows that the unexplained part drop quickly by adding PC1 to PC6, and gradually slower after that. Therefore, using around 6 to 9 PCs would explain the majority of the variance in our data.

To further developed the correlations, we also draw a biplot. According to the principle of the biplots, the variables are positively correlated if their vectors form an angle smaller than 90 degrees. From our plot, all the assets are having a acute angle, so we conclude that they all having a positively relationship.

df = return[,2:13]
model_pca = prcomp(df, center = TRUE, scale. = TRUE)
su = summary(model_pca)
su
## Importance of components:
##                           PC1    PC2     PC3     PC4    PC5     PC6     PC7
## Standard deviation     2.0761 1.3586 1.06307 0.93727 0.8879 0.80884 0.78736
## Proportion of Variance 0.3592 0.1538 0.09418 0.07321 0.0657 0.05452 0.05166
## Cumulative Proportion  0.3592 0.5130 0.60717 0.68038 0.7461 0.80060 0.85227
##                            PC8     PC9    PC10    PC11    PC12
## Standard deviation     0.70899 0.67221 0.56977 0.54539 0.44294
## Proportion of Variance 0.04189 0.03766 0.02705 0.02479 0.01635
## Cumulative Proportion  0.89415 0.93181 0.95886 0.98365 1.00000
screeplot(model_pca, type = "line", main = "Scree plot")

autoplot(model_pca,  loadings = TRUE, loadings.colour = 'blue',
         loadings.label = TRUE)

#model_pca$rotation[,1:2]

factor analysis

X = scale(df,center=T,scale=T)
n=fa.parallel(X)

## Parallel analysis suggests that the number of factors =  2  and the number of components =  2
df.fa = fa(X,nfactors = 4,fm="mle",rotate = "varimax")
fa.diagram(df.fa)

faa1 = factanal(df,factor=3,rotation = "varimax")
faa1
## 
## Call:
## factanal(x = df, factors = 3, rotation = "varimax")
## 
## Uniquenesses:
##  MSFT  TSLA  AAPL  TWTR  AMZN    FB  NFLX   AAL   DAL   BAC  NVDA   WBD 
## 0.459 0.744 0.534 0.853 0.411 0.617 0.626 0.230 0.193 0.502 0.553 0.338 
## 
## Loadings:
##      Factor1 Factor2 Factor3
## MSFT  0.673   0.208   0.210 
## TSLA  0.462           0.190 
## AAPL  0.639   0.201   0.130 
## TWTR  0.148           0.348 
## AMZN  0.752           0.150 
## FB    0.568           0.237 
## NFLX  0.598           0.122 
## AAL   0.148   0.855   0.130 
## DAL           0.880   0.171 
## BAC   0.298   0.488   0.414 
## NVDA  0.634   0.201         
## WBD   0.101   0.220   0.776 
## 
##                Factor1 Factor2 Factor3
## SS loadings      2.870   1.934   1.135
## Proportion Var   0.239   0.161   0.095
## Cumulative Var   0.239   0.400   0.495
## 
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 40.08 on 33 degrees of freedom.
## The p-value is 0.185

6.Risk Management

Return <- return
Price <- asset
n = nrow(Return)
alpha = 0.05
# Non-paramatric

VaR_non = c()
ES_non = c()
for (i in 2:13) {
  q = as.numeric(quantile(unlist(Return[,i]), alpha)) #-0.06918179
  VaR_nonp = -100000 * q
  IEVaR = (Return[,i] < q)
  ES_nonp = -100000 * sum(Return[,i] * IEVaR) / sum(IEVaR)
  options(digits = 5)
  VaR_non = c(VaR_non, VaR_nonp)
  ES_non = c(ES_non, ES_nonp)
}
VaR_non
##  [1]  6918.2 15147.2 10749.9 22137.6  8881.9 10561.0 13072.7 15706.1 12891.9
## [10] 13041.2 11974.8 14903.9
ES_non
##  [1]  9183 20688 13639 25981 12875 15968 20587 24467 20512 16550 21517 19432
# Parametric - normal dist

# m = mean(Return$MSFT)
# sd = sd(Return$MSFT)
# qalpha = qnorm(alpha)
# VaR_par = -100000 * (m + sd * qalpha)
# ES = 100000 * (-m + sd * (dnorm(qalpha))/alpha)

VaR = c()
ES = c()
for (i in 2:13) {
  m = mean(unlist(Return[,i]))
  sd = sd(unlist(Return[,i]))
  qalpha = qnorm(alpha)
  VaR_par = -100000 * (m + sd * qalpha)
  ES_par = 100000 * (-m + sd * (dnorm(qalpha))/alpha)
  VaR = c(VaR, VaR_par)
  ES = c(ES, ES_par)
}
VaR
##  [1]  7239.2 23648.3  9963.7 23401.8 10831.8 11724.8 16743.7 18793.0 14803.4
## [10] 12075.0 14227.4 18533.6
ES
##  [1]  9712 30876 13200 29475 14243 15119 21650 23596 18777 15498 19149 23269

Bootstrap

library(bootstrap)
## 
## Attaching package: 'bootstrap'
## The following object is masked from 'package:pls':
## 
##     crossval
b=1000; n=400; S0=100000; var.boot=rep(0,b); es.boot=rep(0,b)

#non-parametric bootstrapping VaR
CI_nonp = c()
ES_nonp = c()
for (k in 2:13){
  for (i in 1:b){
    data.boot=sample(unlist(Return[,k]),n, replace=TRUE)
    var.boot[i]= -S0 * q
    #ES
    q = as.numeric(quantile(data.boot, 0.05))
    IEVaR = (data.boot < q)
    es.boot[i]  = -100000 * sum(data.boot * IEVaR) / sum(IEVaR)
    options(digits = 5)
    }
    CI_nonp = rbind(CI_nonp, c(quantile(var.boot, 0.025), quantile(var.boot, 0.975)))
    ES_nonp = rbind(ES_nonp, c(quantile(es.boot, 0.025), quantile(es.boot, 0.975)))
}

CI_nonp
##          2.5%   97.5%
##  [1,]  6610.1  7642.1
##  [2,] 14710.9 20336.6
##  [3,]  8108.5 11679.8
##  [4,] 19101.8 26024.1
##  [5,]  7861.3 10283.0
##  [6,]  9596.2 11187.7
##  [7,] 12068.5 19338.2
##  [8,] 14919.7 19849.6
##  [9,] 10980.4 15905.7
## [10,] 11456.9 13239.5
## [11,]  8707.3 18231.9
## [12,] 13540.7 18046.4
ES_nonp
##          2.5% 97.5%
##  [1,]  8062.5 10881
##  [2,] 18381.3 22243
##  [3,] 12208.0 15841
##  [4,] 24251.5 27276
##  [5,] 10885.7 16104
##  [6,] 12315.7 22335
##  [7,] 17424.7 24924
##  [8,] 20042.6 30785
##  [9,] 16429.8 27491
## [10,] 14336.8 20078
## [11,] 17571.3 24602
## [12,] 17529.6 21865
#parametric
CI_p = c()
ES_p = c()
for (k in 2:13){
  for (i in 1:b){
    data.boot=sample(unlist(Return[,k]),n, replace=TRUE)
    var.boot[i]= -S0*(mean(data.boot)+qnorm(0.05)*sd(data.boot))
    #ES
    s = sd(data.boot)*sqrt((1000-1)/1000)
    m = mean(data.boot)
    es.boot[i] = 100000 * (-m + s * (dnorm(qalpha))/alpha)
    }
    sd = sd(var.boot)*sqrt((1000-1)/1000)
    mean = mean(var.boot)
    CI_p = rbind(CI_p, c(mean + sd*qnorm(0.025), mean + sd*qnorm(0.975)))
    ES_p = rbind(ES_p, c(quantile(es.boot, 0.025), quantile(es.boot, 0.975)))
}
CI_p
##          [,1]    [,2]
##  [1,]  6279.0  8100.5
##  [2,] 21236.0 25752.1
##  [3,]  8682.5 11105.3
##  [4,] 21286.4 25179.3
##  [5,]  9567.7 11902.3
##  [6,]  9846.9 13421.5
##  [7,] 14874.3 18447.2
##  [8,] 16834.6 20502.2
##  [9,] 12629.6 16745.7
## [10,] 10566.4 13464.4
## [11,] 12295.1 16018.5
## [12,] 17044.9 19856.0
ES_p
##          2.5% 97.5%
##  [1,]  8574.4 10711
##  [2,] 27918.9 33628
##  [3,] 11648.3 14480
##  [4,] 27096.8 31459
##  [5,] 12741.4 15430
##  [6,] 12972.6 17073
##  [7,] 19407.2 23767
##  [8,] 21294.1 25589
##  [9,] 16312.3 21192
## [10,] 13873.4 17146
## [11,] 16796.7 21127
## [12,] 21350.5 24862

7.Copulas

library(copula)

#12 assets
est.df = data.frame()
i = 2
while (i <= 13) {
  est = as.numeric( fitdistr(unlist(Return[,i]),"t")$estimate )
  est[2] = est[2] * sqrt( est[3] / (est[3]-2) )
  est.row = c(est[1], est[2], est[3])
  est.df = rbind(est.df, est.row)
   i = i + 1
}
colnames(est.df) = c('m', 'sd', 'df')
# est.df: dataframe for estimations
est.df
##              m       sd      df
## 1   0.02427852 0.059806  6.1502
## 2   0.02107033 0.179692  3.8043
## 3   0.02794368 0.077174 85.9747
## 4   0.00336264 0.144514 26.7937
## 5   0.02330902 0.081512  8.4090
## 6   0.01693417 0.080778  5.0603
## 7   0.02101915 0.118101  5.4579
## 8   0.00098099 0.115981 10.4759
## 9   0.00847228 0.094574  4.9511
## 10  0.01587257 0.081778  7.2435
## 11  0.05033777 0.117408  9.9174
## 12 -0.01123768 0.119099  3.9940
data1 = cbind(rank(Return$MSFT)/(n+1), rank(Return$TSLA)/(n+1), rank(Return$AAPL)/(n+1), rank(Return$TWTR)/(n+1), rank(Return$AMZN)/(n+1), rank(Return$FB)/(n+1), rank(Return$NFLX)/(n+1), rank(Return$AAL)/(n+1), rank(Return$DAL)/(n+1), rank(Return$BAC)/(n+1), rank(Return$NVDA)/(n+1), rank(Return$WBD)/(n+1))

#fit t copula
cop_t_dim12 = tCopula(dim = 12, dispstr = "un", df=n-1)
ft1 = fitCopula(data = data1, copula = cop_t_dim12, method="ml")
ft1@estimate
##  [1]  0.85475  0.91235  0.80902  0.90781  0.88447  0.87286  0.85582  0.84566
##  [9]  0.89127  0.89719  0.87230  0.87739  0.82898  0.86659  0.86716  0.86022
## [17]  0.81634  0.80362  0.84942  0.85955  0.84101  0.83140  0.89427  0.89283
## [25]  0.85483  0.85002  0.84778  0.85846  0.89800  0.86533  0.84412  0.84790
## [33]  0.83700  0.79633  0.80647  0.82372  0.81253  0.82629  0.89189  0.91085
## [41]  0.83263  0.82139  0.86491  0.89658  0.84579  0.88087  0.82927  0.81514
## [49]  0.87028  0.85706  0.84999  0.81481  0.79458  0.84251  0.87361  0.83042
## [57]  0.94989  0.89428  0.84669  0.85586  0.88559  0.83657  0.85692  0.84774
## [65]  0.89854  0.82264 52.47837
ft1@loglik
## [1] 1445.5
AIC(ft1)
## [1] -2757.1
# fit gaussian copula
fnorm = fitCopula(data = data1, copula = normalCopula(dim=12), method="ml")
fnorm@estimate
## [1] 0.85056
fnorm@loglik
## [1] 1345.8
AIC(fnorm)
## [1] -2689.6
# frank
ffrank = fitCopula(copula = frankCopula(3, dim = 12), data = data1, method="ml")
ffrank@estimate
## [1] 15.113
ffrank@loglik
## [1] 1395.4
AIC(ffrank)
## [1] -2788.7
# clayton
fclayton = fitCopula(copula = claytonCopula(3, dim = 12), data = data1, method="ml")
fclayton@estimate
## [1] 1.0806
fclayton@loglik
## [1] 1078.3
AIC(fclayton)
## [1] -2154.7
# gumbel
fgumbel = fitCopula(copula = gumbelCopula(3, dim = 12), data = data1, method="ml")
fgumbel@estimate
## [1] 3.4771
fgumbel@loglik
## [1] 1301
AIC(fgumbel)
## [1] -2600.1
# joe
fjoe = fitCopula(copula = joeCopula(2, dim = 12), data = data1, method="ml")
fjoe@estimate
## [1] 13.161
fjoe@loglik
## [1] 1420
AIC(fjoe)
## [1] -2838.1